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CUBIC  SPLINE  SOLUTION  OF  INTEGRAL  EQUATIONS 


I.  INTRODUCTION 


A general  form  in  which  a linear  integral  equation  can  be  written 


is 


fd(x)  g'(x)  + f^x)  8(x)  + fin(x) 


/• 


dy  K(x,y)  g(P(x,y))  , (1 


where  the  functions  f^Cx),  f^^(x),  f^(x)  and  K(x,y)  are  known  and  g(x) 
is  to  be  solved  for.  The  subscripts  d,  li  and  in  stand  for  derivative, 
linear  and  inhomogenous , respectively.  Most  commonly,  such  integral 
equations  are  not  susceptable  to  exact  solution;  some  approximate  or 
numerical  method  must  be  used.  This  report  presents  a version  of  a 
simple  method  based  on  numerical  quadrature  and  spline  interpolation. 
Essentially,  the  spline  interpolation  method  is  used  to  functionally 
represent  g(x)  in  terms  of  a set  of  parameters.  Given  this  parameter- 
ization, we  can  perform  the  various  operations  indicated  in  Eq.  (1), 
obtaining  a simple  matrix  equation  that  can  be  solved  for  the  unknown 
parameters. 

1-4 

The  method  of  spline  interpolation  of  a function  g(x),  in  its 
usual  form,  involves  breaking  an  Interval  [a,b]  into  M-l  subintervals 
and  representing  the  function  g(x)  on  each  subinterval  by  a polynomial 


(we  use  cubic  splines  exclusively).  The  M points  {x^}  that  define  the 


subintervals  are  called  knots;  the  M values  of  g(x)  at  the  knots  are 


customarily  known  (g^gfx^)).  In  the  case  of  cubic  splines,  there  are 


4 (M-l)  parameters  to  be  set— -M  of  them  by  the  values  g^  themselves. 
In  addition,  the  values  of  the  splines  and  their  first  and  second 
derivatives  are  matched  at  the  M-2  interior  knots.  We  are  then  left 


Note:  Manuscript  submitted  December  29,  1978. 
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with  4 (M-l )- [M+3  (M-2)  ] “2  conditions  still  needed  in  order  to  completely 
specify  the  system.  There  are  a variety  of  ways  of  choosing  these  last 
two  condtions;  several  are  discussed  in  Appendix  A. 

In  the  case  of  an  integral  equation,  the  values  g^  are  not  known, 
but  we  can  formally  carry  out  the  spline  analysis  anyway,  obtaining 
splines  that  represent  g(x)  in  terms  of  the  M unknown  values  g^.  We 
then  have  a functional  (piecewise  polynomial)  representation  of  g(x)  and 
can  proceed  to  evaluate  the  integral  in  the  integral  equation  by  an 
appropriate  method  of  quadrature— generally  a numerical  method.  What  is 
left  is  a matrix  equation  of  the  form 


U g - V , 
i.J  8j  i ’ 


where  the  g are  the  unknown  quantities,  U are  the  elements  of  a 
J t J 

square  matrix  and  the  are  the  inhomogeneous  terms  from  Eq.  (1). 

(We  use  the  notation  that  repeated  indices  are  summed  over.)  Thus, 

solving  the  original  integral  equation  is  reduced  to  straightforward 

matrix  operations. 

Methods  similar  to  this  have  been  described  in  the  literature^  ^ 

and  tests  and  qualifications  on  their  use  have  been  investigated.^’^ 

Unfortunately,  we  find  these  results  too  particular  to  the  specific 

cases  investigated  to  predict  from  them  whether  there  are  instabilities 

in  our  use  of  the  method.  Instead,  we  adopt  a more  pragmatic  approach. 

Given  some  particular  integral  equation  (1)  to  be  solved,  we  construct 

a similar  integral  equation  whose  exact  solution  we  know.  Based  on 

intuitive  arguments,  one  often  has  an  idea  of  the  gross  behavior  of  the 

solution  g(x)  of  a particular  integral  equation.  We  choose  a simple 

function  g(x)  that  exhibits  that  gross  behavior.  Further,  we  choose  a 

form  for  K(x,y)  that  approximates  the  real  kernel  but  is  simple  enough 

that  we  can  perform  the  integration  in  Eq.  (1)  exactly.  Similarly,  we 

choose  functions  f^  and  f^  that  approximate  those  in  the  original 

integral  equation.  We  then  set  f so  that  the  given  trial  g(x)  is  the 

in 

exact  solution  of  the  constructed  integral  equation. 

This  stratagem  yields  an  integral  equation  that  can  be  used  to  test 


2 


Che  various  numerical  procedures  chac  we  want  Co  use  on  Che  original 
Integral  equation.  Since  we  tailored  the  properties  of  the  exactly 
solvable  equation  to  the  original,  we  are  reasonably  confident  that 
numerical  procedures  that  work  well  for  the  one  will  not  give  spurious 
values  for  the  other. 

In  the  next  sections  we  describe  our  method  in  detail.  Included  is 
a generalization  of  the  usual  spline  approach  that  allows  known  discon- 
tinuities in  the  derivatives  of  the  solution  g(x)  to  be  built  into  the 
spline  equations.  Following  that  we  give  an  example  of  the  sort  of 
model  equations  that  we  solve  in  order  to  test  our  procedures.  For 
pedagogical  purposes,  we  present  in  Appendix  B a simple,  exactly  solved, 
example  of  our  method. 


II.  INTEGRAL  EQUATION  AND  NOTATION 


The  type  of  integral  equation  that  we  will  solve  using  a cubic 
spline  representation  of  the  solution  and  numerical  quadrature  to  per- 
form the  necessary  integrations  is 


fd(x)  g'(x)  + fj^Cx)  g(x)  + fin(x) 


2(x) 

(x) 


dy  KCx.y)  g(P(x,y) ) , 


(3 


which  is  similar  to  Eq.  (1).  The  various  functions  are  defined  on  the 
interval  [a,b] . The  limits  of  integration. 


a<y1(x)<y2(x)<b  , (4 

could  be  incorporated  into  the  kernel  K(x,y),  but  are  more  convenient  in 
this  explicitly  stated  form.  Similarly,  we  presume  that  P(x,y)  is  a 
continuous  function  of  x and  y such  that 


a<P(x,y)<b  . (5 

For  the  moment,  we  further  presume  that  functions  f^fx),  f^(x) 


and  f (x),  as  well  as  P(x,y)  and  the  kernel  K(x,y),  are  continuous 
in 

with  continuous  first  and  second  derivatives. 


III.  CUBIC  SPLINE  INTERPOLATION 
A.  The  Spline  Equations 

Before  concentrating  on  the  problem  of  solving  Integral  equations, 
let  us  review  the  method  of  cubic  spline  interpolation.  Our  notation 
follows  that  of  Ref.  4.  Suppose  that  som*  function  g(x)  is  known  at  M 
points  on  the  Interval  [a,b] , 

a»x1<x2<. • •<xM_^<xM“b  . (6 

These  M points  are  known  as  knots.  In  each  interval  [x^.x^j],  we 
represent  g(x)  by 

h 2 1 

8i(x)  " 8i  9i(x)  + »i+i  “ 61  fgi  {si  “ si  > 

+ g^!  (t£  - tt3)]  . 1-1, M-l  (7 

where 

hi  ■ x14-1  - , 1-1 , M-l , (8 


st  - s^x)  - (*1+1  - x)/ht  , 1-1, M-l  , 


(9 


C1  “ “ (x  - Xj^/hj^  , 1-1, M-l. 


(10 


The  M values  g^g^,  * • • ,g^j  are  the  second  derivatives  of  g(x)  evaluated 
at  the  knots  x^  ; they  are  determined  by  the  continuity  conditions  at 
the  knots. 

In  particular,  the  form  (7)  was  chosen  so  that 


8<xl)  “ gt  “ 


(11 


4 


(12 


g"<x1)  - g'^  - g'^x^  - g'^U^  ; 

that  Is,  the  values  and  second  derivatives  of  adjacent  cubic  functions 
are  equal  at  each  knot.  By  setting 

s;<*i)  * * i”! .M— 1 (13 

we  obtain  the  further  condition 


Vi  «;.i + 2(hi-i + v + hi  ®i+i 

“ & f(g1+1  “ gi^hi  “ “ 8i-l^hi-l^  » 1*2, M-l  (14 

which  gives  M-2  equations  for  determining  the  M unknown  values  g^'  in 
terms  of  the  g^.  Two  more  conditions  are  needed;  for  the  moment  we 
will  use 


*1  " *2  3nd  Ci  " *M  * (15 

A discussion  of  these  and  similar  conditions  is  given  in  Appendix  A. 

Henceforth,  we  will  use  the  notation  that  each  index  runs  from  l 
to  M,  unless  otherwise  indicated.  Also,  repeated  indices  are  understood 
to  be  summed  over,  from  1 to  M.  We  can  express  Eq.  (14)  in  matrix  form 


i"  - 


J.i  1 


6 V *i 


where  all  elements  of  L vanish  except 


(16 


Ll,l  “ 

Hl,M  " 1 

Ll,2- 

Hi, M-i  " _1 

Li,i  " 

2(hi-l  + V 

L . . 

■ h,  , 

1,1-1 

i-1 

L 

■ h 

•*4,1+1 

i 

1-2.M-1 


(18 


Similarly,  all  elements  of  the  matrix  R vanish  except  for 


Ri,i  “ -1/hi-l  -1/hi 
Ri,i-1  " 1/hi-l 


i-2,M-l  . 


(19 


i,i+l 


1/h, 


Upon  multiplying  Eq.  (16)  by  the  inverse  of  the  matrix  L,  we  obtain 


.-1 


*;  ■ 6 a h.j  Rj,k 


(20 


We  define  two  new  matrices 


Ai,k  " _hi  (L  }i,j  Rj,k  » 1"l,M“1 


Bi,k  “ “hi  (L  ^i+l.j  Rj,k  * i’"1»M“1 


(21 

(22 


so  that  we  can  write  the  spline  functions  in  the  form 


8i(x)  " 8i  3i  + «1+l  ci  + Ai,k  «k  (si  ’ 8i  > 

+ Bi,k  gk  (tt  “ tt  ) . i-l.M-1. 


(23 


In  solving  the  integro-dif f erential  equation,  we  shall  also  need 
the  derivatives  of  g(x)  at  the  knots.  By  taking  the  derivative  of  each 
spline  function  we  obtain 

8i(xi)  " [ 8i+l  8i  + 2 Ai,k  8k  + Bi,k  8k  ]/hi  ’ 1“1,M_l  (24 

8M-1(xM)  [ " ^-1  “ *11-1, k 8k  " 2 ®M-l,k  8k  1 /hM-l  * (25 

which  could  be  rewritten  in  the  form 


8i  ‘ °i,k  8k  ’ i“1*M-1- 


(26 


6 


B.  Inclusion  of  Discontinuities  in  Derivatives  of  g(x) 

The  solutions  of  the  integro-dif f erential  equations  in  which  we  are 
interested  often  have  discontinuities  in  their  first  or  second  deriv- 
atives. The  cubic  spline  method  is  readily  adapted  to  account  for  such 
discontinuities. 

When  we  create  a mesh  of  knots  for  the  cubic  spline  procedure,  we 
naturally  set  a knot  at  each  of  the  points  of  discontinuity.  For  nota- 
tional  simplicity,  we  will  suppose  that  the  function  g(x)  possesses 
discontinuities  in  its  first  and  second  derivatives  at  each  of  its 
knots.  We  write  these  as 


Gi  " yio  dx  (g(x+y)  - 8(x-y))lx,x  , i-2,M-l  (27 

and  similarly  for  the  second  derivative  G'^.  For  the  moment,  these 
discontinuities  are  presumed  to  be  known.  The  addition  of  these  discon- 
tinuities generates  the  following  changes  in  our  equations:  The  contin- 
uity equations  for  derivatives  (12)  and  (13)  become 

g£_l(Xi)  + G^(Xl)  - , i»2,M-l  (28 

g^_l(Xi)  + G'^)  - g^Xj)  » i“2 ,M-1 . (29 

The  constants,  such  as  g^,  that  appear  in  the  expressions  for  the  spline 
functions  are  the  values  of  derivatives  of  g(x)  as  determined  from  below 
the  knots;  for  example,  g^-g^^).  Thus  Eq.  (7)  for  g^x)  becomes 

h2  3 > 

gf(x)  - gt  st(x)  + g1+1  t^x)  - -£i  [(g^  + G£)  (st  - s^) 

+ (ti  - <3° 


and  Eq.  (14),  which  arises  from  the  continuity  equation  for  the  first 
derivatives,  becomes 


7 


* **  S * 


* * 

\ 

I'  * 

jv' 


r 


Mm  - 


(31 


*,• 

< 


hl-l  gi-l  + 2 (hl-l  + hi)  *i  + hl  «l+1  ' -2  hi  G'l  ‘hl-l  Gi-1 


-6  G 


[ + 6 ((g1+1  “ 8i)/hi  " , 1-2.M-1 


The  definitions  (17-19)  of  the  matrices  L and  R do  not  change,  but 
Eqs.  (16)  and  (20)  become 


L a"  » 6 R b + r 

J.i  8i  J.i  gi  i 


-1 


(L  ’i.J  [ 6Rj,k  gk  + rJ]  * 


(32 

(33 


where 


ri  - -6  G'  - 2 h±  G'{  - hJL_1  G'^_1  . 


(34 


If  auxiliary  conditions  other  than  (15)  are  used,  slight  changes  in  the 
definitions  of  both  the  matrices  L and  R and  the  vector  r will  be  needed 
(see  Appendix  A).  With  the  use  of  the  matrices  A and  B,  defined  in 
Eqs.  (21)  and  (22),  and  with 


Xt  - - -|  hi2  ( (L”1)1  k rk  + G^  ] , i-l.M-l 


Yi  “ “ 6 hi  (L  ^i+l.k  rk  * i”1,M“1 


(35 

(36 


we  have  our  final  expression  for  the  spline  functions. 


gi(x)  - gt  st  + g1+1  tt  + (Sl  - 3l  ) [Ai>k  gk  + 


+ (ti  “ V > [Bi,k  gk  + V • 


<*7 


With  the  addition  of  the  discontinuities,  the  first  derivative  of 
the  spline  functions  evaluated  at  (just  above)  each  knot  can  be  written 


8I(xi)  “ °i,k  gk  + di  * i“1*M’1 


(38 


i 


k. 


where 

°i.k  = 

[ 5k,i+l  " 5k,i 

°M,k 

[ &M,k  " &k,M-l 

di  * [ 

2 X + Yt  ] / h 

dM“  [ 

l 

3? 

1 

1 

ro 

3 

1 

(39 

(40 

(41 

(42 


IV.  SOLVING  THE  INTEGRO-DIFFERENTIAL  EQUATION 
A.  Evaluation  of  the  Integral  Term 

Let  us  rewrite  our  basic  integral  equation,  Eq.  (3),  with  the  vari- 
able x taking  on  only  the  values  of  the  knots  x^.  We  then  have 


fd(xt)  g(xt)  + fIt(x,)  g(x1)  + 


(43 


where  for  each  of  the  M knots  there  is  a different  integral 


fy 2(xi) 

Jy  i<*i> 


dy  K(xt,y)  g(P(x1>y)) 


(44 


to  be  evaluated.  Further,  we  choose  for  each  value  of  x^  a numerical 
quadrature  procedure,  with  associated  weights  w^  ^ and  ordinates  y^  k, 
so  the  1^  can  be  expressed  as 


(45 


We  suppose,  in  order  to  keep  the  notation  simple,  that  each  quadrature 
uses  N points;  actually,  an  entirely  different  quadrature  procedure 
could  be  used  for  each  x^. 

For  each  value  of  P(X£»y^  we  determine  the  quantity  p(i,k) 
such  that 


9 


• *m-  ■*'  A- 


(46 


Xp(i,k)  ± P(xi^i,k)  ± 


p(i,k)+l  ’ 


so  that  we  can  replace  g(P(x1>yi  ) by  the  appropriate  spline 

functKm  » >k)(p(VTl>k)).  / 


We  rewrite  1^  in  the  form 


I.  =*  I.  g + J.  , 
i i,m  °m  i * 


(47 


where 


Ji  ' ”i,k  K<VJri,k)  1 XP  (sp  - %3> + YP  (tP  - 'p3>  1 • <48 

h.j  ’ "i,k  K(V’'i,k)  1 *P,j  *P  * Vuj  CP 

+ , (s^  - s 3)  + B , (t  - t 3)  ] . (49 

P.j  P P P.j  P P 

For  given  i and  k,  p**p(i,k)  is  determined  by  Eq.  (46);  s and  t are 

P P 

defined  by 


3P  ■ (xp+i  - p(xi^i,k))  ' hp  * 


cP  = (p(xi’yi,k)  - y ' \ • 


(50 

(51 


B.  Matrix  Form  of  the  Integral  Equation 

We  can  now  collect  the  terms  generated  in  the  last  few  sections  and 
rewrite  the  integral  equation  in  the  form 


Ui,m  «m  “ Vi  * 


which  has  the  solution 


gm  - (iT1)  , V . 

m m,  i i 


(52 


(53 


10 


i/. 


* 

i 


The  quantities  U and  V in  Eq.  (52)  are 


I,m  I,m  11  i l,m  d 1 


Vi  " - Ji  + fin(xi)  + fd(xi}  di  * 


(54 

(55 


V.  TEST  CASE 


As  a numerical  test  of  our  method,  we  use  an  integral  equation 

that  models  a radiation  damage  process.  Specifically,  as  a heavy  ion 

traverses  an  amorphous  material,  some  of  its  energy  is  lost  to  elastic, 

displacement  producing,  collisions.  The  following  integral  equation 

yields  the  energy  g(x)  that  the  ion  with  initial  energy  x gives  up  to 

8“  1 0 

such  elastic  collsions.  The  integral  equation  has  the  same  form 

as  Eq.  (3);  specifically. 


fd(x)  g'(x)  + f-^x)  8(x)  + ftn(x) 


g(y)  U(y-l) 


. h -k 

ay  x * y * 


g(x-y)  U(y-l) 


(56 


with  the  add^d  condition  that  g(0)«0.  The  various  functions  in  Eq.  (56) 
are  given  in  the  following  table;  U(y)  is  the  unit  Heavyside  function. 


Functional  Forms  Appearing  in  Equation  (56) 


1 

Function 

0<x<l  | 

l<x<2 

I 2<x 

1 

1 

1 

fd  1 

1 

2 \e 

kx  + 2x^  | 

kx2+2x^ 

| kx2+2x  2 

1 

1 

1 

fli  1 

o 1 

2(x^  -1) 

1 2(xs  -1) 

1 

1 

f in  1 

-kx2  - 2x^  | 

fL(x) 

| f1(x)+f2(x) 

1 

11 


Z 


r*» 


- 


where  and  f?  are  given  by 

, 3 , 3 , 

f^(x)“-2kx  +0.5kx  2-x^lnx-4x2  +10(x-x^  )+l  (57 

f2(x)«  2x 2 [2(x+l)+(4-5x) (x-1)  ^ -sin  ^x  2 +sin  *((x-l)/x)^  ] (58 

The  exact  solution  of  this  integral  equation  is  g(x)-x  for  0<x<l 
and  g(x)“2x-x2  for  x>l. 

With  the  condition  that  g(0)“0,  we  can  examine  Eq.  (54)  for  discon- 
tinuities, discovering  that 


G'(l)«l/2 

G"(l)-l/4  (59 

and  that  there  are  no  other  discontinuities.  The  computer  program  that 

we  wrote  to  solve  this  sort  of  integral  equation  has  a grid  of  48  values 

of  x , ranging  from  x-0  to  a somewhat  arbitrarily  chosen  x-x  “400. 
i max 

The  integrations  are  done  using  Guass  quadrature,  with  a change  of  vari- 
able to  compensate  for  the  y 2 behavior  of  the  integrand.  (More  par- 
ticularly, the  interval  from  0 to  x_>v  is  broken  into  a number  of  sub- 
max 

intervals,  with  different  quadratures  used  in  each.) 

The  result  of  this  procedure  is  that  the  spline  solution  of  the 

integral  equation  reproduces  the  exact  results  to  no  worse  than  0.3%-- 

generally  the  agreement  is  much  better.  Of  course,  the  error  increases 

and  decreases  as  x is  increased  and  decreased.  The  quality  of  this 
max 

agreement  leads  us  to  be  confident  in  using  the  same  methods  f^r  similar 
integral  equations  whose  exact  solutions  are  not  known. 


APPENDIX  A:  AUXILIARY  CONDITIONS 


As  we  noted  in  the  introduction,  the  straightforward  application  of 
cubic  spline  interpolation  requires  the  setting  of  two  auxiliary  condi- 
tions to  be  satisfied  by  the  spline  coefficients.  Usually,  one  of  these 


' 

< 


is  imposed  at  each  end  of  the  interval,  so  that  the  various  matrices 
remain  balanced  about  their  diagonals.  Below,  we  list  five  of  the 
conditions  that  we  have  used  at  x-x^,  noting  that  for  each  there  is  a 
corresponding  condition  that  can  be  used  at  x-x^. 

We  use  the  notation  ^ * [ ] to  indicate  the  components  of 

the  first  row  of  the  H matrix;  any  elements  in  that  row  that  are  not 
listed  are  taken  to  be  zero.  Similar  notation  is  used  to  indicate  the 
elements  of  a vector. 

The  conditions  are: 

a)  Suppose  that  the  second  derivative  of  g(x)  is  known  at  x^, 
so  that  g"(x^)»g*,  then  the  auxiliary  condition  is 

Ll,i  " C1’J  ’ Rl,i  “ KM*  ri  • (A1 

b)  If  the  first  derivative  is  known  at  x ^ , say  g'(x^)«g*,  then  by 
differentiating  Eq.  (7)  we  obtain 

6 (g2  ■ gj)/h1  * 6 g*  ■ 2 hj  g”  + hj  g£  (A2 

LUl  - [2ht  ,hx , ] , - t-l/h1,i/h1>],  v±  “ [-6g* , ] . (A3 

c)  If  both  g'(x^)  and  g"(x^)  are  known,  we  can  use  (a)  and  (b), 
provided  that  Eq.  (30)  is  not  used  for  i-2  and  provided  that  G£  and  GlJ 
are  included,  if  necessary. 

d)  In  the  absence  of  any  particular  information,  we  can  set 

gj  “ a*  g’^  , (A4 

Ll,i  " » Rl,i  “ [0‘]»  T±  *10*1  • (as 

where  a*  is  some  constant.  We  usually  take  a*«l,  but  1/2  and  0 are  also 
commonly  used.^ 

e)  We  can  also  set  the  third  derivatives  of  g(x)  equal  at  x-Xj, 
so  that 
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(A6 


h2  8"  " (h1+h2>  82  + g'j  - 0 , 

i “ [hj ,h j “ [0,],  r^  “[0,]  . (A7 


APPENDIX  B:  A SIMPLE  EXAMPLE 


In  this  section  we  present  a simple  example  of  our  method,  by  solv- 
ing an  integro-dif ferential  equation  whose  exact  solution  we  know.  We 
will  outline  the  solution  by  supplying  the  various  matrices  and  vectors. 
We  start  with  the  equation 

-(x-2)  g * (x)  + (x2-3x+8)  g (x)  + (-2x3+6x2-14x-4) 

- 3 (x-2) J 2*  dy  U(y-2)  g(y)  (B 1 

where  U(y)  is  the  unit  step  (Heavyside)  function.  We  choose  the  knots 
xi  ■ 1,2, 3, 4 , so  that  h -1,  1-1,3. 

If  we  examine  Eq.  (Bl),  and,  in  particular,  differentiate  it  twice, 
we  find  that 

g(2")  - g<2+)  - 4 

g'(2“)  - g'(2+)  - 2 
g"(2*)  - 0 

g"(2+)  - 6,  (B2 

so  that  all  of  the  G'  and  G"  quantities  are  zero  except 

G"2  - 6.  (B3 

This  example  is  simple  enough  that  all  of  the  evaluations  can  be 
done  exactly,  including  the  integrals  over  the  spline  polynomials 
in  Eq.  (Bl).  We  thus  can  successively  form  various  the  matrices  and 
vectors  that  are  defined  in  Sections  III  and  IV.  They  are: 
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> 

V 

t 


l 

t 

f 


| 1 -1  0 

o 1 

1 o 

0 

0 

o 1 

114  1 

o 1 

1 1 

-2 

1 

o 1 

L * 1 0 1 4 

1 1 

R - 

1 o 

1 

-2 

1 1 

|0  0-1 

1 1 

1 o 

0 

0 

o 1 

(B4 

r - [0,  -12,  -6,  0] 

| 19  5 

-1  1 1 

-1  1 -5  5 
24  L " | 1 -1 

-1  1 1 

5 -5  | 

1 1 -1 

5 19  | 

(B5 

| 5 -11 

7 -1  1 

1 5 

-11 

7 

-1  1 

I 5 -11 

7 -1  1 

1 -1 

7 

-11 

5 1 

“24  A ' | -1  7 - 

•11  5 | 

-24 

B - 

1 -1 

7 

-11 

5 1 

1 “ 

— — | 

1 — 

— 

— 

— | 

(B6 

8 X - [3, -5,1,—] 

8 

Y - 

[3,1,1 

, — ] 

(B7 

1 -13  19-7  1 | 

1 o 

0 

0 

o 1 

I -3  -3 

7 -1  1 

1 o 

0 

0 

0 I 

8 D " | 1-7 

3 3 | 

16 

I - 

1 -2 

26 

26 

-2  I 

1 -1  7 -19  13  | 

1 -2 

38 

122 

34  | 

(B8 

8 d - [9,-9, 3,-3]  16  J - [0,0, -6, -6]  (B9  •, 

The  vector  forms  of  the  various  functions  In  Eq.  (Bl)  are 
fd  - [1,0, -1,-2]  fu  - [6,6,8,12]  fln  - [-14,-24,-46,-92]  (BIO 


so  that  the  vector  V is  * 


The  matrix  U and  its  inverse  are 


| -70  -38  14  -2  | 

| 0 -96  0 0 | 

16  U “ | 0 12  -96  4 | 

| -6  66  46  -106  | 


(B12 


16374  U 


-3747 

1453 

-522 

51  1 

0 

-2729 

0 

o 1 

9 

-423 

-2778 

-105  | 

(Bl  3 

216 

-1965 

-1176 

-2520  1 

Upon  multiplying  V by  U 


we  obtain 


gt  “ [2,4,9,20] 


(B14 


If  we  use  these  values  of  in  Eq.  (37)  for  the  spline  func- 
tions, we  obtain 

g^x)  - 2x 

g2(x)  - 2x  + 3(x-2)2  (B 15 

g3(x)  - g2(x)  . 

Thus  we  exactly  reproduce  the  function 

g (x ) - 2x  + 3(x-2) 2 U(x-2)  (B16 

that  we  used  to  create  Eq.  (Bl).  Of  course,  the  reason  the  solution  is 
exact  is  that  only  simple  polynomials  are  involved,  but  we  see  the  ease 
with  which  discontinuities  can  be  handled. 
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